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We describe a numerical method to compute Casimir forces in arbitrary geometries, for arbitrary 
dielectric and metallic materials, with arbitrary accuracy (given sufficient computational resources). 
Our approach, based on well-established integration of the mean stress tensor evaluated via the 
fluctuation-dissipation theorem, is designed to directly exploit fast methods developed for classical 
computational electromagnet ism, since it only involves repeated evaluation of the Green's func- 
tion for imaginary frequencies ( equivalent ly, real frequencies in imaginary time). We develop the 
approach by systematically examining various formulations of Casimir forces from the previous 
decades and evaluating them according to their suitability for numerical computation. We illustrate 
our approach with a simple finite-difference frequency-domain implementation, test it for known 
geometries such as a cylinder and a plate, and apply it to new geometries. In particular, we show 
that a piston- like geometry of two squares sliding between metal walls, in both two and three di- 
mensions with both perfect and realistic metallic materials, exhibits a surprising non-monotonic 
"lateral" force from the walls. 



I. INTRODUCTION 

One of the most dramatic manifestations of quantum 
mechanics observed in the last half-century is the Gasimir 
force: a tiny force on an uncharged, source-free body due 
to changes in the zero-point energy associated with quan- 
tum vacuum fluctuations (virtual photons) [H [2] [3l [H |5] . 
There have been many experimental verifications of the 
Gasimir force reported in recent decades [H [71 E] , but al- 
ways restricted to simple geometries (parallel plates [9], 
spheres and plates [a[IOl[IIl[Il[T3[Tl[l3[ni[IIl[IHl 
El EQl EB 122], or crossed cylinders [lOl [m |23] ) . More- 
over, the force in these particular geometries is almost 
always attractive [24 (except possibly for some unusual 
material systems [i|2ll[2i[27l[2il23[30l|3ll[3^ 
and monotonically decreasing with separation. Thus, one 
might ask whether it is possible to obtain non-monotonic 
or even repulsive forces in more complex structures, and 
more generally whether complex geometries might give 
rise to unexpected force phenomena. For more com- 
plicated geometries, however, calculations become ex- 
tremely cumbersome and often require drastic approx- 
imations, a limitation that has hampered experimental 
and theoretical work beyond the standard geometries. 

In this paper, we explore several ways in which well- 
established, efficient techniques from standard compu- 
tational electromagnetism can be brought to bear on 
this problem, in order to predict forces for arbitrary ge- 
ometries and materials with arbitrary accuracy (no un- 
controlled approximations). Starting from the simplest, 
most direct approaches, we show that practical consid- 
erations naturally lead towards a particular method in- 
volving the integral of the Minkowski stress tensor by 
repeated evaluation of the imaginary-frequency Green's 
function — a method previously developed for purely an- 
alytical calculations [3, 34, 35 . We illustrate the method 
by a simple finite-difference implementation, but evalu- 
ation of the Green's function is so standard that many 



more sophisticated techniques are immediately applica- 
ble, and we discuss what techniques are likely to be 
optimal. Our approach is tested for geometries with 
known solutions, and then is applied to new geometries 
in two and three dimensions that lead to surprising non- 
monotonic effects. We also demonstrate the applica- 
tion of our technique to dispersive dielectric materials, 
not just for idealized perfect metals. We explain how 
the same technique can be used for visualization of the 
Gasimir interactions between bodies, as well as for com- 
puting other quantities of interest, such as torques. The 
key advantage of exploiting standard computational ap- 
proaches is not merely that existing code, error analy- 
ses, and other experience can be applied to the Gasimir 
problem, but also that these methods have been proven 
to scale to large three-dimensional problems, which have 
previously seemed out of reach of exact methods for the 
Gasimir force. 

The most common approach to predicting the Gasimir 
force has been to consider approximations for small 
perturbations around known solutions, such as parallel 
plates or dilute gases. For parallel plates in d dimen- 
sions, separated by a distance a, there is a well-known 
attractive force that scales as l/a^+^, first predicted by 
Gasimir [1 and later extended to formulas for any planar- 
multilayer dielectric distribution e{x^u;) via the gener- 
alized Lifshitz formula [36^. A direct, intuitive exten- 
sion of this result is the proximity force approximation 
(PFA) [37], which treats the force between two surfaces 
as a two-body interaction given by the sum of "paral- 
lel plate" contributions. Valid in the limit of small cur- 
vature, PFA provides an easy way to conceptualize the 
Gasimir force in complex geometries as a simple two-body 
force law, but unfortunately it may also be deceptive: 
outside its range of applicability, the Gasimir force is 
not additive [2 and may be qualitatively different from 
PFA's predictions H IMl HHl SIl 1121 HSl |44]. Other 
perturbative approaches include renormalized Gasimir- 
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Polder [45l [46] or semi-classical interactions [47], mul- 
tiple scattering expansions ^48| |49|, classical ray optics 
approximations [50], higher- order PFA corrections [51], 
and other perturbative techniques [57. The ray-optics 
approach is especially interesting because, although it is 
only strictly valid in the small-curvature limit, it captures 
multiple-body interactions and can therefore sometimes 
predict the qualitative behavior in cases where other ap- 
proximations fail [SSIISIIES]- Nevertheless, none of these 
methods provide any guarantees of accuracy in arbitrary 
geometries, where they involve uncontrolled approxima- 
tions. Therefore, for complex new geometries, where one 
might hope to encounter behaviors very different from 
those in the parallel-plate limit, a different approach is 
required. 

To this end, researchers have sought "exact" numeri- 
cal methods applicable to arbitrary geometries — that is, 
methods that converge to the exact result with arbi- 
trary accuracy given sufficient computational resources. 
One such method was proposed by Ref. [56l based on a 
path-integral representation for the effective action; this 
method has predicted the force between a cylinder and 
a plate [57 , and between corrugated surfaces [40l [41]. 
It is based on a surface parameterization of the fields 
coupled via vacuum Green's functions, requiring 0{N'^) 
storage and 0{N^) time for N degrees of freedom, mak- 
ing scaling to three dimensions problematic. Another 
exact method is the "world-line approach" [58], based 
on Monte-Carlo path-integral calculations. The scaling 
of the world-line method involves a statistical analysis, 
determined by the relative feature sizes in the geome- 
try, and is discussed below. The methods of Ref. 56 
and Ref. 58 have currently only been demonstrated for 
perfect-metallic z-invariant structures — in this case, the 
vector unknowns can be decomposed into TE (E • z = 0) 
and TM (H • z = 0) scalar fields with Neumann and 
Dirichlet boundary conditions, respectively — although 
generalizations have been proposed [58l |59]. Here, we 
propose a method based on evaluation of the mean 
Minkowski stress tensor via the fluctuation-dissipation 
theorem, which only involves repeated evaluation of the 
electromagnetic imaginary-frequency Green's function. 
For our initial volume discretization with N degrees of 
freedom and an efficient iterative solver, this requires 
0{N) storage and at best 0{N'^~^^^) time in d dimen- 
sions. Furthermore, because evaluation of the Green's 
function is such a standard problem in classical com- 
putational electromagnetism, it will be possible to ex- 
ploit many developments in fast solvers, based on finite- 
element [HOl [nH [Ml [63], spectral [60l [64], or boundary- 
element methods [60l [61] [62] [65] [66] [67]. As we argue 
below, a future implementation using boundary-element 
methods should attain nearly 0(7VlogA^) time. To il- 
lustrate the method, however, our initial implementation 
is based on the much simpler finite-difference frequency- 
domain method [68] with a conjugate-gradient solver [69] . 
as described below. 

In the following sections, we describe the step-by-step 



conceptual development of our computational method. 
Our purpose here is to start back at the beginning, with 
the earliest theoretical descriptions of the Casimir force, 
and analyze these formulations from the point of view 
of their suitability for purely numerical calculations. Al- 
though the final technique, in terms of the stress ten- 
sor integrated over space and imaginary frequency, can 
be viewed as a numerical implementation of a textbook 
result due to Dzyaloshinskii et al. [S] [35l [70], it is il- 
lustrative to derive it as the culmination of a sequence 
of simpler approaches, in order to show how it circum- 
vents a number of numerical obstacles that would hin- 
der a more direct method. To begin with, we illustrate 
the methods using the well-known case of parallel plates 
where they can be compared to analytical expressions, 
but a more rigorous test is subsequently provided by the 
situation of a cylinder and plate, recently solved numeri- 
cally [57]. Finally, we apply our method to a new geom- 
etry of a "piston" -like structure involving blocks sliding 
between parallel walls, in both two and three dimensions 
with both perfect metals and more realistic dispersive 
dielectrics, and demonstrate a surprising non-monotonic 
"lateral" effect from the walls. We conclude by analyz- 
ing the scaling of the method compared to previous ap- 
proaches and discussing the application of more sophis- 
ticated finite-element and boundary-element techniques. 



II. A SIMPLISTIC APPROACH 

Perhaps the simplest approach to computing the 
Casimir force is to think of it as the derivative of the 
zero-point energy U expressed as a sum of ground-state 
photon energies hjo/2. For each photon mode with fre- 
quency cj, the zero-point energy is fujo/2^ and thus the 
total Casimir energy, at least for non-dissipative systems 
where uj is real [36] , is formally given by the sum over all 
modes |1][2]: 

U = Y,\hw (1) 

This sum is formally infinite, because the classical har- 
monic modes uo have unbounded frequencies. There is 
some controversy over the physical interpretation of this 
divergence ^1], but in practice it is removed by regular- 
izing the sum in some fashion, for example multiplying 
by e~*^ for 5 > 0, and taking s ^ only after the sum is 
differentiated to obtain the force F = —dU/da between 
two bodies with separation a [72^. This approach, which 
was an early method to analytically compute the force 
between perfect-metal plates [T, might seem to provide 
the most direct computational method as well. After all, 
the computation of electromagnetic eigenmodes is rou- 
tine even in very complicated geometries, and efficient 
methods are known [60l [61] [62l [73]. Unfortunately, it 
turns out not to be practical for this problem (except 
in one-dimensional geometries [74|), as explained below. 
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but the reason why it is impractical points the way to 
more efficient methods. 

To ihustrate the difficulty in directly evaluating 
Eq. ([T]), let us consider the simplest one-dimensional ge- 
ometry: two parallel perfect-metal plates, separated by 
a distance a, in which case one can predict analytically 
an attractive force F = T:hc/2Ao? [71 . Ignoring this 
analytical result, let us apply a numerical method that, 
conceptually, we could apply to an arbitrary geometry: 

1. Discretize space with resolution Ax using a finite- 
difference approximation, with space truncated to a 
finite computational cell (e.g. with periodic bound- 



can be found analytically: 
2 



- 



(3) 



for n = 0, ...,(i/Ax. The energy U is then given by 
summing cOn in Eq. ([T]) for d = a and d = L — and the 
force F by the discrete derivative of U as above. 

Applying this procedure numerically for Ax = 0.05a, 
one obtains the correct force to within 5% (and to any 
desired accuracy by increasing L and decreasing Ax), so 
at first glance it may seem that the method is successful. 
However, its impracticality is revealed if we examine the 
contribution of each frequency Un to the force. 



2. Solve numerically for the eigenmode frequencies co 
and sum to obtain U{a). 

3. Shift one body (one plate) by one pixel Ax and 
thus compute U{a-\- Ax). 

4. Obtain the force F ^ -[U{a + Ax) - /7(a)] /Ax. 

Note that this method automatically provides its own 
regularization: the number of modes in a discretized 
computational cell is finite (the frequencies are bounded 
by the Nyquist limit), and hence U is finite for Ax > 
0. The periodic boundaries will lead to artificial "wrap- 
around" forces, but since Casimir forces decay rapidly 
with distance, the contribution of these forces can be 
made negligible for a sufficiently large computational cell. 
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FIG. 1: (Color) Schematic of one dimensional geometry, show- 
ing two Id metal plates separated by a distance a, embedded 
in a computational cell of length L = a + 4a, with periodic 
boundary conditions, and resolution Ax. 

This method, for the one-dimensional parallel-plate ge- 
ometry, is illustrated in Fig. [l] Here, we have two plates 
with separation a and an overall computational cell size 
L = 5a (which will contribute an erroneous wrap-around 
force 1/16 of the physical force). Maxwell's equations, in 
one dimension, can be written as the scalar eigenproblem 



(in c = 1 units), which is discretized to 



E, 



n+l 



2En + Ern 



Ax2 



^^En 



(2) 



in a center-difference approximation for Ez{nAx) = En. 
For two metal plates with separation d and Dirichlet 
boundary condition Ez = 0^ the discrete eigenvalues cOn 
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FIG. 2: (Color) top: Plot of force summand, or spectral den- 
sity, hALj/2Ax vs. 00 for Id parallel metal plates from Fig.[l] 
bottom: Plot of force partial sum hAuj/2Ax vs. u. 



The top panel of Fig. |2] shows the contribution 
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hAuj/2Ax [Au = uj{a -\- Ax) — uj{a)] to the force sum- 
mation, and the bottom panel shows the corresponding 
partial sum (for frequencies up to co). We see that every 
frequency (of the regularized/finite-resolution problem) 
makes a non-negligible contribution to the force, and the 
summation is of a wildly oscillating quantity that leaves 
a tiny remainder at the end. The reason for these oscilla- 
tions is quite simple: as a is increased, the frequencies on 
the 4a side of the plates increase slowly, while the smaller 
number of frequencies on the a side of the plate decrease 
more rapidly, and these lead to the positive and negative 
contributions in the top panel of Fig. [2j respectively. 

These two features, which are intrinsic properties 
not limited to this particular discretized geometry [75] . 
combine to make this method impractical in higher- 
dimensional structures. Because every frequency con- 
tributes to the force, in a numerical method one must 
compute every eigenvalue of the Maxwell eigenproblem. 
In one dimension, that is not so bad, but in general if 
there are N degrees of freedom {N grid points), then 
computing every eigenvalue of an x matrix requires 
0{N'^) storage and 0{N^) time. This is impractical in 
three dimensions where N may be in the millions. Fur- 
thermore, the wild oscillations of the summand imply 
that the eigenvalues must be computed quite accurately, 
and may exacerbate numerical difficulties in larger prob- 
lems. 

However, these undesirable features are avoidable, be- 
cause we have not yet exploited a key property of 
Maxwell's equations: causality. If we ignore the causality 
constraint, then the oscillatory spectrum would be an ob- 
servable effect: one would simply employ a material that 
is a perfect metal in some frequency range and trans- 
parent otherwise, in order to obtain the force spectrum 



where Gjk solves: (V x V x — cc;^£)G/c(cc;; x — x') = 
(5^(x — x')e/c, with e/e denoting the unit vector in the kth 
direction. For non-dissipative systems in which Eq. ([T]) is 
valid, e is real and we can therefore pull the Im outside of 
the integral. (The generalization to dissipative materials 
is discussed below.) 



Furthermore, we know from causality requirements 
that the Green's function has no poles in the upper half 
plane in cj-space [2j [76]. This means that one can per- 
form a contour integration to relate duo to the integral 
dw along the imaginary- frequency axis co = iw^ also 
known as a Wick rotation. We therefore obtain: 



integrated only in that range. Such a material, how- 
ever, would violate the Kramers-Kronig constraints that 
follow from causality considerations Thus, we are 

motivated to exploit causality in some fashion to avoid 
the oscillatory spectrum. 



III. WICK ROTATION AND ENERGY 
DENSITY 

In order to exploit causality, we will rewrite Eq. ([T]) in 
terms of the electromagnetic Green's function via an inte- 
gral over the density of states. Causality implies that the 
Green's function has no poles in the upper-half plane, so 
one can perform a contour integration, or Wick rotation, 
to transform the sum over real frequencies into an integral 
along the imaginary-frequency axis. The result of this 
standard trick turns out to be a well-known expression: 
an integral of the mean electromagnetic energy density, 
evaluated by the fluctuation-dissipation theorem using 
the temperature (Matsubara) Green's function. Again, 
we focus on the method's suitability as a purely numer- 
ical approach, for arbitrary geometries, and we will find 
that it still leaves something to be desired. 

First, we can express the zero-point energy of Eq. ([T]) 
as an integral over the local density of states I)(x,a;): 

/7 = ^ ^ dco J u;D{^, co) d^x. (4) 

Since we are solving for the eigenstates of Maxwell's equa- 
tions (V X V X —uj'^e)'E = 0, the local density of states 
I)(x, cj) can be expressed in terms of the Green's tensor 



(5) 

I 



U=^ dw — — ^trG(z^;x-x)d^x, (6) 

2n Jo J dw ^ ' ^ ' V ; 

where the new problem to be solved is that of finding the 
solutions to the imaginary-time Green's function (c = 1 
units): 

[V X V X ^w'^eiiw, x)] Gk{iw] x - x') = ^^(x - x')e/e. 

(7) 

As usual, this is formally infinite, because the Green's 
function is singular at x = x', but one typically regu- 
larizes the problem by subtracting the vacuum Green's 
function, which removes the singularity without changing 
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k=l 



V X V X 



|x;ek) = 



1 dju^e) 

TT dcj 



Imtr G(cj; x — x) 



5 



the net force. 

Equation (|6| is not a new result, nor is it limited to 
non-dissipative materials (unlike our derivation) [3l [36]. 
In fact, it is equivalent to the mean energy in the fluctu- 
ating electromagnetic field, derived from the fluctuation- 
dissipation theorem via the temperature Green's func- 
tions [3 . Our purpose in deriving it this way is to empha- 
size the connection to the simplistic approach of Eq. ([T]) . 
In particular, the mean energy in the electromagnetic 
fields (for the case of non-magnetic materials /i = 1) is 
given by [H [76] : 



U = 



dw — 
J 2v 



dw 
d{w'^e) 
2w dw 



(8) 



using the fact that J £(E^)^ = J(H^)^. Here, the key 
point is that the mean values of the fluctuating fields are 
given, via the fluctuation-dissipation theorem, in terms 
of the imaginary- frequency Green's function [3]: 



(i?,(x)i?fe(x')). 



h 



w^Gjkiiw] X — x'). 



(9) 



where the dyadic Green's function Gij solves Eq. ^ 
and obeys the usual boundary conditions on the electric 
field from classical electromagnetism Substituting 
Eq. ([9| into Eq. (|8|, one recovers Eq. |6|. 

From a computational perspective, the imaginary- 
frequency integral of Eq. (|6| turns out to be greatly su- 
perior to the real- frequency summation of Eq. ([T]), for 
two reasons. First, while every real frequency a; con- 
tributed to the force —dU/da^ the same is not true for 
the derivative of the imaginary-frequency integrand. In 
particular, as discussed below, the force integrand in 
imaginary frequencies turns out to be a smooth, non- 
oscillatory, strongly peaked function of w^ meaning that 
one can integrate it via a smooth-quadrature method that 
evaluates the integrand at only a small number of w val- 
ues. Second, the imaginary-frequency Green's function 
turns out to be quite easy to obtain by relatively stan- 
dard methods, including for dissipative systems (where 
obtaining the eigenmodes is harder because it involves a 
non-Hermitian eigenproblem) . These two favorable fea- 
tures are closely related. 

Consider Eq. ^ for the imaginary-frequency Green's 
function. Unlike its real- frequency counterpart, the lin- 
ear operator on the left-hand side of this equation is 
real- symmetric and positive- definite (for w > {)). This is 
true even for dissipative materials where e^uo) is complex, 
since causality requirements imply that £{iw) is purely 
real and positive (in the absence of gain) [76]. For 
one thing, this implies that the most powerful numeri- 
cal methods are applicable to solving the linear system 
of Eq. ^ — many of these methods (e.g. the conjugate- 
gradient method) are restricted to Hermitian positive- 
definite operators ^69^. Also, the resulting Green's func- 



tion is particularly well-behaved: it is exponentially de- 
caying and non-oscillatory. This transforms the highly 
oscillatory real-cj force integrand into a mostly non- 
oscillatory integrand, and also makes the force integrand 
exponentially decaying for large w (for large w^ the in- 
teractions between bodies become exponentially small). 
(In addition, as we will discuss in Sec.|Vin" the exponen- 



tially decaying Green's function is especially favorable for 
boundary-element numerical methods.) 

Again, considering the simplest possible finite- 
difference scheme, this leads us to the following numerical 
algorithm to compute the force: 

1. For a given imaginary frequency w. 

(a) For every grid point x, solve Eq. ([t]) for each 
polarization k to obtain Gkkij^w] x — x). 

(b) Sum over x to compute the spatial integral in 
Eq. (|6|. 

(c) Repeat the above for the body shifted by one 
pixel Ax and subtract to obtain the force in- 
tegrand at w. 

2. Employ a smooth quadrature scheme to integrate 
the above function over w to obtain the force. 

Again, the spatial discretization provides its own regular- 
ization {Gkk is finite), and thus no additional regulariza- 
tion is required (the contribution of the vacuum Green's 
function to the net force is zero). Again, one can truncate 
the computational cell in a number of ways, for example 
with periodic boundaries, and the artifacts thereby intro- 
duced will decrease rapidly with cell size. Again, there 
are also many other ways that one could potentially solve 
for Gkk besides a finite-difference approximation, but we 
will delay discussion of those techniques until we have 
formulated the final method in the next section. 

The above procedure can again be applied to the one- 
dimensional problem of the force between two plates, 
as in Fig. [2| to illustrate its basic features. The re- 
sult is shown in Fig. |3) plotting the force integrand as 
a function of imaginary frequency w^ and the difference 
from Fig. [2] is striking. The integrated force is the same 
(correct) result as before. Unlike Fig. [2j the force inte- 
grand has no sign oscillations, is exponentially decaying 
for large w^ and is strongly peaked around a character- 
istic w = 27rc/a, corresponding to a "wavelength" of a 
(the separation) . These features imply that the force can 
be accurately integrated by an adaptive Gauss-Kronrod 
quadrature scheme [78 using at most a few dozen w 
points. 

Although this method is much more efficient than the 
one described in the previous section, and is potentially 
practical at least in two dimensions, it still has some un- 
desirable features. Suppose that we have N grid points 
in our discretized operator {N may be very large in 
3d). Even if we have an ideal iterative solver for the 
sparse linear system of Eq. ^ , such as an ideal multigrid 
solver [63, l79|, each evaluation of the Green's function 
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FIG. 3: (Color) Plot of Casimir force integrand dU/da be- 
tween two Id parallel plates separated by a distance a = 1 
versus imaginary frequency w = Imcj, using the method of 
Sec. lull 



takes at best 0{N) time with 0{N) storage. However, 
we must evaluate the Green's function 3N times in or- 
der to perform the spatial integration, resulting in 0{N'^) 
complexity. As is discussed in the next section, one can 
do much better than this by using the stress tensor in- 
stead of the energy density. In fact, as is discussed in 
Sec. VIII[ it should ultimately be possible to obtain the 
force with nearly 0(A^log7V) work using a boundary- 
element method to compute the stress tensor, in which 
is only the number of degrees of freedom required to 
represent the interfaces between materials. 



IV. STRESS-TENSOR COMPUTATIONAL 
APPROACH 

After having analyzed the feasibility of several tech- 
niques to solve the Casimir problem through the lens of 



numerical electromagnetism, we are ready to appreciate 
and explore the most feasible of the methods thus far pre- 
sented: an approach based on the Maxwell stress tensor. 



As derived by Dzyaloshinskii et al. [3| |35l [70] , the net 
Casimir force on a body can be expressed as an integral 
over any closed surface around the body of the mean elec- 
tromagnetic stress tensor (T^j), integrated over all fre- 
quencies. Again, using the same arguments as above, it 
is computationally convenient to perform a Wick rota- 
tion, expressing the net force as an integral over imagi- 
nary frequencies uo = iw. (The original derivation used 
imaginary frequencies to start with, via the tempera- 
ture Green's function, but the result is equivalent to a 
Wick rotation of the real-frequency expression. A re- 
lated analytical treatment, but using purely real uo and 
therefore unsuitable for numerical computation because 
of the oscillations discussed above, has also been exam- 
ined [^[81].) The resulting net force is: 



poo 

dw (U) (T(r,i^)) -dS. 

^ surface 



(10) 



In two or one dimensions, one or two of the spatial in- 
tegrals are omitted, respectively, but the result still has 
the units of force; the change in dimensions of dS is bal- 
anced by a change in the dimensions of the delta func- 
tion in Eq. On the other hand, for a 3d structure 
that is z-invariant [a constant 2d cross-section e{x^y)] 
or ?/z-invariant [a univariate e{x)]^ the integrals over the 
invariant directions are replaced by integrals over the cor- 
responding wavevectors, resulting in a net force per unit 
length or per unit area, respectively. (These wavevector 
integrals are discussed in more detail in Sec. VI ) The 
stress tensor is defined as usual by: 



(11) 



where ja and e are the relative permeability and permit- 
tivity, respectively, although in most cases we set /i = 1 
for simplicity (since most materials have negligible mag- 
netic response at short wavelengths, and in any case the 



stress tensor is normally evaluated over a surface lying in 
vacuum). As before, the connection to quantum mechan- 
ics arises from the correlation functions of the fluctuat- 
ing fields, given via the fluctuation-dissipation theorem in 
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terms of the imaginary-cj dyadic Green's function. The 
correlation function for the electric field (EiEj) is given in 
Eq. In this case, however, we also need the magnetic- 
field correlation functions, which can be obtained by dif- 
ferentiating the electric-field Green's function 

{H,{r)Hj{r')) = --(Vx),,(Vx),^G,^(z^; r - ^ , 

TT 

(12) 

(Alternatively, the {HiHj) correlation function can be 
computed from the magnetic Green's function, which is 
the magnetic field in response to a given magnetic-dipole 
current.) The above expressions are given at zero temper- 
ature; the nonzero-temperature force is found by chang- 
ing J dw in Eq. (10) into a discrete summation [3l [49j . 
Although the Green's function (and thus T) is formally 
infinite at r = r', this divergence is conventionally re- 
moved by subtracting the vacuum Green's function; in 
a numerical method with discretized space, as below, 
there is no divergence and no additional regularization 
is required. (The vacuum Green's function gives zero net 
contribution to the dS integral, and therefore need not 
be removed as long as the integrand is finite.) 



Historically, this stress-tensor expression was used to 
derive the standard Lifshitz formula for parallel plates, 
where Gij is known analytically. Its adaptability and 
suitability as a purely computational method does not 
seem to have been recognized, however. As in the previ- 
ous section, the method involves computing the Green's 
function for many imaginary frequencies w and spatial 
points X, integrated over w and x. However, a quick 
glance at Eq. (10) will suggest at least two obvious com- 



putational advantages compared to the method discussed 
in Sec. |III| First, in framing the problem in terms of 
the stress tensor, we have reduced the spatial integral 
over the whole volume (Eq. ([8|) to a surface integral 
around the body of interest. This implies that, for a 
(i-dimensional geometry, the computational effort due to 
spatial integration is reduced from 0{N'^) to 0{N'^~^^^). 
Second, the force is now given directly in terms of the 
dyadic Green's function (via the stress tensor), rather 
than its derivative, which avoids another layer of compu- 
tation. Moreover, although our derivation is only valid 
when the stress tensor is evaluated at points within loss- 
less dielectrics (regardless of whether the bodies them- 
selves are dissipative) , one can also extend it for evalua- 
tion in absorbing media [70 . However, the case discussed 
above (bodies separated by vacuum) is the most common. 



So far, we have presented the step-by-step development 
of an efficient approach to computing Casimir forces. 
In what follows, we illustrate our new approach using 
a proof-of-concept finite-difference implementation, and 
present results that demonstrate its fiexibility and utility. 



V. THE FINITE-DIFFERENCE METHOD 

At this point, all that remains is the numerical com- 
putation of the Green's function via Eq. ^ for an imag- 
inary frequency co = iw. This is simply the inversion of a 
linear operator [V x V x -\-w'^£{t, iw)] that has the con- 
venient properties of being real-symmetric and positive- 
definite, as stated above. Almost any technique devel- 
oped for computational electromagnetism is applicable 
here, modified to operate at an imaginary frequency. 
To illustrate our approach, we used a very simple, yet 
extremely general, method: finite-difference frequency- 
domain (FDFD) discretization of Eq. ([t]) in a staggered 
Yee grid [68], which we then invert by a conjugate- 
gradient method [69]. Although the Yee grid in princi- 
ple allows second-order-accurate finite-difference approx- 
imations, unfortunately the whole scheme becomes only 
first-order-accurate once a discontinuous dielectric func- 
tion e is included. (There are ways to treat interfaces 
more accurately [82], but we did not implement them 
here.) Moreover, a very high resolution is often required 
to resolve the stress tensor close to a dielectric boundary 
due to the Green's function divergences as a boundary is 
approached [65 . Despite its shortcomings, however, we 
found FDFD to be sufficient to obtain accurate results 
(to within a few percent in a reasonable time) for two- 
dimensional, and three-dimensional z-invariant, geome- 
tries. The key advantage of FDFD is its ffexibility: with 
very little effort, we were able to implement support for 
arbitrary geometric shapes and arbitrary materials (both 
perfect metals and dispersive/dissipative dielectrics). 

Again choosing the simplest possible approach, we ap- 
ply periodic boundary conditions at the edges of the 
computational cell, which are accurate as long as the 
boundaries are sufficiently far compared to the separa- 
tion between the interacting bodies. That is, the period- 
icity leads to artificial "wrap-around" forces that decay 
rapidly with cell size L (at least as in 2d); we chose 
cell sizes large enough to make these contributions neg- 
ligible (< 1%). 





FIG. 4: (Color) Schematic illustration of a possible contour 
around a body; the force on the body is given by an integral 
of the stress tensor around this contour. 

The computational process (using a simple finite- 
difference scheme) goes as follows: 

1. Pick a contour /surface around the body of interest, 
as in Fig. [4] (which will typically not coincide with 
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the boundary of the body). 
2. For a given frequency w: 

(a) For every grid point x on the discretized con- 
tour/surface, solve Eq. ^ for each polariza- 
tion k to obtain Gjk{iw; x — x). 

(b) Integrate the resulting stress tensor Tjk over 
the surface, as in Eq. (10). 



3. Integrate the above function over w to obtain the 
force; since the integrand is a smooth function of 
an efficient adaptive quadrature scheme can be 
employed [78 . 

Although this scheme does not require any additional 
regularization (the integrand is finite for a finite resolu- 
tion, and the integral of the vacuum stress tensor over 
the contour is zero), we have found that numerical con- 
vergence can be accelerated by subtracting the stress- 
tensor integral over the isolated bodies. For example, 
in the schematic of Fig. |4| we would first compute the 
stress-tensor integral for the two bodies as shown, then 
subtract the integral of the stress tensor over the same 
surface with one body removed, and then subtract again 
for the stress tensor with the other body removed. Of 
course, these subtracted quantities are zero in the limit 
of infinite spatial resolution — there is no net force on an 
isolated body. However, at a finite resolution the dis- 
cretization error at the interface between two materials 
can lead to a finite force that vanishes as resolution is 
increased. By subtracting this error term from the force, 
we find that the numerical error is typically reduced by 
an order of magnitude or so. We emphasize, however, 
that this is merely an optimization — even without sub- 
traction, the force converges to the correct result, and 
merely requires a somewhat higher resolution. 



VI. FORCES IN TRANSLATION-INVARIANT 
STRUCTURES 

It is common to solve for the Casimir force between 
bodies that are translation-invariant in one or more di- 
rections; for example, between a cylinder and a plate [57] 
that are invariant in the z direction. More generally, one 
might consider structures that are periodic in some direc- 
tion with a non-zero period A, where A ^ corresponds 
to translation invariance. Intuitively, in such cases one 
need only perform computations in the unit cell of the 
periodicity, and the spatial integration in the invariant 
direction (s) is replaced by integration over a wavevector 
k from Bloch's theorem [83l|84]. Although special cases 
of this familiar idea are well known in Casimir computa- 
tions [56 , 84 , here we provide a review of this approach 
for an arbitrary periodicity in the context of stress-tensor 
computational methods; the detailed derivation is pro- 
vided in the Appendix. We also mention a useful opti- 
mization for the common special case of perfect-metallic 
z-invariant structures. 



Let us consider a single direction of periodicity: sup- 
pose that the structure is periodic in z with period A. In 
this case, it is natural to choose a surface for our stress- 
tensor integral that is also periodic in z. For example, 
imagine that Fig. [4] depicts a two-dimensional {xy) cross- 
section of a z-invariant structure, and the dashed line 
depicts a cross-section of the corresponding 2;-invariant 
stress-tensor surface. Because the total force is infinite, 
the quantity of interest is the force per unit z. It is con- 
venient to consider the net force F from a finite length 
L = NA with periodic boundaries, and obtain the force 
per unit length as limAr^ooF/L. Naively, F/L can be 



written directly via Eq. (10), where we break the inte- 



gral over z into a summation over the unit cells: 



F _ 1 

L~ L 

_ 1 

~ A 



/ V rrT(z^;r-nAz) -dS, (13) 

•^0 n=0 



dw 



JJ T(i^;r) • dS, 



(14) 



where the surface integral is over the portion of the sur- 
face lying in the unit cell only, and in the second line 
we have used the fact that the stress tensor is periodic. 
This expression is inconvenient, however, because the di- 
rect evaluation of T{iw^ r) requires the response to a sin- 
gle point source in the large-L structure, and a single 
point source does not produce a periodic field (requiring 
a full three-dimensional calculation even for a z-invariant 
structure). Rather, one would like to consider the field 
in response to periodic point sources, which produce a 
periodic field that can be treated by a small computa- 
tional cell with periodic boundary conditions. This is 
accomplished by Fourier-transforming the expressions in 
Eq. ( 13 ) and taking the N ^ oo limit, as described in de- 



tail by the Appendix. The resulting force per unit length 



is: 



F 

L 



1 
A 



poo pn/A 77 



where the surface integral is still over the portion of the 
surface lying in the unit cell. Here, T{iw,kz;r) denotes 
the stress tensor computed from the Green's functions 
for Bloch-periodic boundaries — that is, from the fields in 
response to a periodic set of point-dipole sources with 
phase e*^^^^ in the nth unit cell. This stress tensor can 
be computed using a computational cell that is only one 
unit cell in the z direction, e.g. by a two-dimensional 
computational cell for a z-invariant structure. [Equiva- 
lently, T{iw^kz]T) could instead be computed from the 
Green's function for ordinary periodic boundaries, but 
with V replaced by V + ikzZ [83 .] Just as for the 
stress tensor is a smooth function of kz and therefore the 
kz integral can be computed by an efficient quadrature 
scheme (e.g. Gaussian quadrature). 

If the structure is periodic (or invariant) in more than 
one direction, one simply repeats the above procedure: 
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for each periodic direction, we only consider the por- 
tion of the stress-tensor integral in the unit cell, with 
Bloch-periodic boundary conditions, and integrate over 
the corresponding Bloch wavevector component. Also, 
by symmetry one only need integrate over the irreducible 
Brillouin zone of the structure [83], e.g. in one dimen- 
sion (where time-reversal symmetry normally equates 

kz and —kz) the Jl^^j^dkz integral can be replaced by 

For the common case of a z-invariant perfect-metal 
structure (i.e. one has a homogeneous e surrounded by 
perfect metal), there are several important simplifica- 
tions. First, the {kz^w) Green's function is exactly the 
same as the (0, y^^^Fh-^) solution [57 . Therefore, we 
need only compute the kz = solutions at each and 
weight the dw integrand by a factor of ttw (the circum- 
ference of a semi-circle of radius w) divided by the 27r 
that would appear in the dkz integral. (For two direc- 
tions of translational symmetry, one would weight the 
integral by the area of a hemisphere of radius and in 
general, in d dimensions, by the area of a radius-i(; hy- 
persphere.) Furthermore, at /c^ = 0, the solutions can be 
divided into two scalar polarizations, TE (E • z = 0) and 
TM (H-z = 0) [76]. 




0.5 1 1.5 2 2.5 3 3.5 4 

CO = iw (271 c/a) 

FIG. 5: (Color) Plot of Casimir force integrand ■ dS be- 
tween two Id parallel plates separated by a distance a = 1 
versus imaginary frequency w = Imcj, using the stress-tensor 
method of Sec. [iVj 



B. Cylinder and plate 



VII. NUMERICAL RESULTS AND 
VISUALIZATION 

In the following sections we demonstrate our method's 
validity by checking it against known results for perfect 
metals, and in particular for the case of a cylinder ad- 
jacent to a plate. A new geometry that displays an in- 
teresting non-monotonic behavior is presented for both 
perfect and realistic dispersive metals, and in both 2d 
and 3d. Furthermore, we explain how one can use stress- 
tensor maps to visualize the interactions between bodies 
and identify the most important spatial regions. 



A. Parallel plates in one dimension 

First, for comparison to the one-dimensional parallel- 
plate integrands plotted in Fig. |2] and Fig. |3j we plot the 
corresponding integrand for the stress tensor in Fig. |5j 
Again, this integrates to the correct force (7rfic/24a^ for 
each polarization), but we note that the integrand is not 
identical to the integrand from the imaginary-frequency 
energy derivative. 

Since the one-dimensional parallel-plate force is com- 
monly derived from the Lifshitz formula, which in turn is 
derived from the stress tensor, this cannot be regarded as 
a rigorous validation of our method (except in the most 
basic sense of checks for bugs in our code). 



A more complicated geometry, consisting of a perfect 
metallic cylinder adjacent to a perfect metallic plate in 
three dimensions, was solved numerically by Ref. i57i to 
which our results are compared in Fig. [6| Ref. ^ used 
a specialized Fourier-Bessel basis specific to this cylin- 
drical geometry, which should have exponential (spec- 
tral) convergence. Our use of a simple uniform grid was 
necessarily much less efficient, especially with the first- 
order accuracy, but was able to match the Ref. |57j re- 
sults within ~ 3% using reasonable computational re- 
sources. A simple grid has the advantage of being very 
general, as illustrated below, but other general bases with 
much greater efficiency are possible using finite-element 
or boundary-element methods; the latter, in particular, 
could use a spectral Fourier basis similar to Ref. 57 and 
exploit a fast- mult ipole method or similar 0{N\ogN) 
solver technique. Surface discretizations (boundary el- 
ements) will also have the advantage that the infinite 
amount of space surrounding the objects is treated an- 
alytically rather than having to be truncated with some 
boundary conditions (here, periodic). This is discussed 
in greater detail in Sec. |VIII[ 

Also shown, in the inset of Fig. |6j is a plot of the 
interaction stress-tensor component A{Txx) at a typical 
imaginary frequency w = 27rc/a. By "interaction" stress 
tensor A(Tjj), we mean the total (Tij) of the full ge- 
ometry minus the sum of the (T^j)'s computed for each 
body in isolation. Here, the stress tensors of the iso- 
lated cylinder and plate have been subtracted, giving us 
a way to visualize the force due to the interaction. As 
described further below, such stress plots reveal the spa- 
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FIG. 6: (Color) Casimir force between a 3d radius- i? cylin- 
der and a plate (inset), relative to the proximity- force ap- 
proximation FpFA, vs. normalized separation a/R. The solid 
lines are the Casimir force computed in Ref. [57lfor TE (gray) 
and TM (blue) polarizations, along with results computed by 
our method with a simple finite- difference discretization (gray 
squares). Error bars were estimated for some data points 
by using computations at multiple spatial resolutions. Inset 
shows interaction stress tensor A{Txx) at a typical imaginary 
frequency w = 27rc/a, where red indicates attractive stress. 



tial regions in which two bodies most strongly affect one 
another, and therefore reveal where a change of the ge- 
ometry would have the most impact. (In contrast, Ref. 85 
plots an interaction-energy density that does not directly 
reveal the force, since the force requires the energy to be 
differentiated with respect to a. For example, Ref. l85l s 
subtracted energy density apparently goes nearly to zero 
as a metallic surface is approached, whereas the stress 
tensor cannot since the stress integration surface is arbi- 
trary.) 



C. Two-dimensional metal piston and 
non-monotonic "lateral" forces 

We now consider a more complicated geometry in 
which there are interactions between multiple bodies: a 
two-dimensional "piston" -like structure, shown in Fig.[7| 
consisting of two metal s x s squares separated by a 
distance a from one another (here, s = a) and sepa- 
rated by a distance h from infinite metal plates on either 
side. We then compute the Casimir force between the 
two squares, in two dimensions (that is, for 2:-invariant 
fields, unlike the cylinder case above where z oscillations 
were included), as a function of the separation h. The re- 
sult for perfect conductors is shown in Fig. [Tj plotted for 
the TE and TM polarizations and also showing the total 
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h/a 

FIG. 7: (Color) Casimir force between 2d (z-invariant fields) 
metal squares F/Fpfa, vs. distance from metal plate h (in- 
set), normalized by the total force (TE+TM) obtained using 
the PEA, FpFA = hcC{3)s/87va^. The total force is plotted 
(black squares) along with the TE (red dots) and TM (blue 
circles) contributions. 



force. (Error bars are not shown because the estimated 
error is < 1%. This structure is computationally easier 
than the cylinder-plate problem, for a finite- difference 
discretization, because the metallic walls allow the com- 
putational cell size to be small in at least one direction.) 
In the limit of — > 0, this structure approaches the 2d 
"Casimir piston," which has been solved analytically for 
the TM polarization [86^. Our results, extrapolated to 
= 0, agree agree with this analytical result to within 
3% (although we have computational difficulties for small 
h due to the high resolution required to resolve a small 
feature in FDFD). For > 0, however, the result is sur- 
prising in at least two ways. First, the total force is 
non-monotonic in /i, due to a competition between the 
TE and TM contributions to the forces. Second, the h 
dependence of the force is a lateral effect of the parallel 
plates on the squares, which would be zero by symmetry 
in PEA or any other two-body-interaction approxima- 
tion. 

The reader may notice that the TE and TM forces in 
the cylinder-plate case. Fig. [6) also have opposite-sign 
slopes in the graph, and one may therefore suspect that 
non-monotonic forces could occur in that case as well. 
However, in the cylinder-plate case this apparent differ- 
ence in sign is merely an artifact of the normalization: 
the PEA force varies with the separation a, and when 
the actual force (which is monotonically decaying with 
a) is divided by this variable force one can obtain "non- 
monotonic" plots. (A similar "non-monotonic" plot can 
be seen in Ref. 57 ) In contrast, for Fig. [7j the PEA nor- 
malization is constant because a is fixed, and thus the 
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relative forces frora different h are directly comparable. 

Although lateral forces can still arise qualitatively 
in various approximations, such as in ray optics or in 
PFA restricted to "line-of-sight" interactions, it may 
not be immediately clear how these could predict non- 
monotonicity. We also note that, in the large- limit, the 
force remains different from PFA due to finite-5 "edge" 
effects [85 , which are captured by our method. It turns 
out that one can qualitatively predict the non-monotonic 
behavior, due to the competition between TE and TM 
forces, using the ray-optics approximation, although this 
approximation is not quantitatively accurate except for 
h = 0;we will describe this ray-optics analysis in a future 
publication [55] , 
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FIG. 8: (Color) (a-f): TM stress map of the geometry in 
Fig. [7| for various h. The interaction stress tensors {Txx) 
(left) and (T,^) (right) for: (a),(d): h = 0.5a; (b),(e): 
h — a] and (c),(f): h — 2a^ where blue/white/red = repul- 
sive / zero /attractive. 



To further explore the source of the /i-dependence, we 
plot the TM interaction-stress maps A (Txx) and A{Txy) 
in Fig. [8} for the perfect-metal squares at a typical fre- 
quency w = 27rc/a, and for varying distances from the 
metal plates {h = 0.5, 1.0, 2.0). As shown, the mag- 
nitudes of both the xx (a-c) and xy (d-f) components 
of the stress tensor change dramatically as the metal 
plates are brought closer to the squares. For example, 
one change in the force integral comes from T^y, which 
for isolated squares has an asymmetric pattern at the 
four corners that will contribute to the attractive force, 
whereas the presence of the plates induces a more sym- 
metric pattern of stresses at the four corners that will 
have nearly zero integral. This results in a decreasing TM 
force with decreasing h as observed in Fig. [Tj Because 
stress maps indicate where bodies interact and with what 



signs, it may be useful in future work to explore whether 
they can be used to design unusual behaviors such as 
non-additive, non-monotonic, or even repulsive forces. 



D. Two-dimensional dielectric pistons 

Our method is also capable, without modification, of 
handling arbitrary dielectric materials. The calculation of 
general dispersive media can be performed with minor 
or no additional computational effort, since the compu- 
tations at different w are independent. Unlike most pre- 
viously published techniques, which do not easily gen- 
eralize to non-perfect metals, the stress tensor approach 
does not distinguish between the two regimes, and com- 
putational methods for inhomogeneous dielectric materi- 
als are widely available. Furthermore, we reiterate that 
along the imaginary-o; axis, e is purely real and positive 
even for dissipative materials (which have complex s on 
the real-cj axis), greatly simplifying computations. 

The method's ability to handle dielectric structures is 
demonstrated below, where the Casimir force between 
the two squares is shown for two different cases: 
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FIG. 9: (Color) Casimir force between 2d (z-invariant fields) 
dielectric (s = 4) squares F/Fpfa, vs. distance from metal 
plate h (inset), normalized by the total force (TE+TM) ob- 
tained using the PFA (Here, the PFA force is computed for 
x-infinite slabs of dielectric £ = 4). The total force is plotted 
(black squares) along with the TE (red dots) and TM (blue 
circles) contributions. 

First, we compute the force between two squares made 
of dielectric material with £ = 4 (an artificial mathemati- 
cal choice for illustration purposes), whereas the parallel 
plates are still perfect metal. The result is shown on 
the plot of Fig. |9] As might be expected, the dielectric 
squares have a weaker interaction than the perfect-metal 
squares, but are still non-monotonic. 
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tor of irw, and therefore the 3d force is obtained from the 
2d calculations with very little additional computation. 




FIG. 10: (Color) Casimir force between 2d (z-invariant fields) 
gold squares F/Fpfa, vs. distance from metal plate h (inset), 
normalized by the total force (TE+TM) obtained using the 
PFA. (Here, the PFA force is computed for x-infinite gold 
slabs). The total force is plotted (black squares) along with 
the TE (red dots) and TM (blue circles) contributions. 



Second, as a more interesting example, the squares are 
made of gold with a Drude dispersion taken from exper- 
iment, again with adjacent perfect metallic plates. In 
particular, the following Drude model is used for the ma- 
terial dispersion of gold [87]: 



CO {u -\- iVp) 



(16) 



with LJn 



1.37 X 10^^ Hz and F 



5.32 X 10^^ Hz, 
corresponding to ujp = 7.2731 and Tp = 0.028243 in our 
units of 27rc/a, for a = 1 /im. For uj = iw^ this is real 
and positive, as expected. The resulting force is shown 
in Fig. [To] Not surprisingly, the gold squares have a 
weaker interaction than perfect-metal squares, since at 
large w = Imco the dielectric constant s goes to 1. 



E. Three-dimensional piston 

The previous 2d calculations are important in at least 
two ways: first, they allow us to check the stress tensor 
method against previous piston calculations in the h ^ 
limit, while exploring an interesting new geometry; sec- 
ond, based on their results, one might predict a similar 
behavior for the force per unit length in the analogous 3d 
2;-invariant geometry, shown in the inset of Fig. 11 In- 
deed, this is the case: the non-monotonic force in three 



dimensions is shown in Fig. [TT] As discussed in Sec. [VTl 
the integrand for a z-invariant perfect-metallic structure 
differs from the two-dimensional integrand only by a fac- 
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FIG. 11: (Color) Casimir force per unit length between z- 
invariant 3d metal blocks F/Fpfa, vs. distance from metal 
plate h (inset), normalized by the total force (TE+TM) ob- 
tained using the PFA, Fpfa = /ics7r^/480a^. The total force 
is plotted (black squares) along with the TE (red dots) and 
TM (blue circles) contributions. 

Again, in the h = limit there are known analyti- 
cal solutions for this geometry based on the ray-optics 
method [88 or the zeta-function technique [89] . Linearly 
extrapolating our plot to = 0, we reproduce these re- 
sults to within ^ 2-3%. 



VIII. BEYOND FINITE-DIFFERENCE 
METHODS 

Above, we implemented the stress-tensor integration 
using a finite-difference frequency-domain approach to 
compute the Green's function. While sufficient for a 
proof-of-concept implementation, one would like to use 
more sophisticated methods in order to explore complex 
geometries more quickly, especially in three dimensions. 
The primary drawbacks of the finite-difference scheme 
are threefold. First, the material discontinuities imply 
that the error converges only linearly with resolution, al- 
though there are techniques to improve this to quadratic 
convergence [82 . Second, non- uniform resolution would 
be desirable to handle small features, such as the nar- 
row channels in the piston structure for small h. Third, 
the stress-tensor integrand is a discretized, non-smooth 
function of space, meaning that we must evaluate it at 
a number of grid points proportional to the resolution 
(or resolution squared, in three dimensions); in general, 
the integrand must therefore be evaluated 0{N^^~^^^^) 
times for a d — 1 dimensional surface in d dimensions. 
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leading to at best 0{N'^~^^^) complexity. To address 
these drawbacks, we consider two standard approaches 
to solving partial differential equations in a more effi- 
cient manner in complex geometries: finite-element and 
boundary-element methods. 

Finite-element methods can employ a non-uniform vol- 
ume discretization, via an unstructured mesh, that can 
both put more resolution where it is needed and conform 
to the material interfaces to obtain higher-order accu- 
racy [60l inU |62l [63]. However, from the perspective of 
Casimir-force calculations, finite-element methods seem 
to have two potential drawbacks. First, because space is 
still discretized, the stress tensor is again not a smooth 
function of space and its accurate integration requires 
that the Green's function be evaluated at many mesh 
points. Second, there may be a problem with regular- 
ization: although the Green's function does not diverge 
in discretized space, with a non-uniform resolution this 
effective regularization varies at different points. Unless 
there is a way to locally regularize the problem (sub- 
tracting a vacuum Green's function computed at the lo- 
cal spatial resolution), this may lead to unphysical, non- 
convergent forces. 

Boundary-element methods (BEMs) involve a dis- 
cretization in terms of unknowns only at the interfaces 
between different materials — these surface unknowns are 
coupled to one another via the (known) Green's func- 
tions of the homogeneous regions [60 , 61, 62 , 65 , 66l [67]. 
They are thus ideal for open problems, in which the bod- 
ies are surrounded by infinite volumes of empty space, 
since those infinite regions are treated analytically. Using 
the fast-multipole method (FMM) or other fast integral- 
equation methods [60l |62j [65] , BEMs can solve for 
the surface unknowns, and hence the Green's function, 
in 0{N\ogN) time for N discretized unknowns, multi- 
plied by a number of iterations (<C N) that depends on 
the condition number of the matrix and the precondi- 
tioning [69 . Furthermore, BEMs have two unique ad- 
vantages when applied to the problem of Casimir forces. 
First, the regularization can be performed analytically: 
since BEMs express the Green's function as the sum 
of the vacuum Green's function plus a set of contribu- 
tions from surface currents, the vacuum Green's function 
can be trivially subtracted analytically. Second, because 
space is not discretized, the stress tensor in a BEM will 
be a smooth (infinitely differentiable) function of space — 
this means that the spatial integral can be performed 
with exponentially-convergent smooth quadrature meth- 
ods. Therefore, the number of times that the Green's 
function must be computed is determined only by the 
convergence of the smooth multidimensional quadrature 
in 2 + 1 dimensions (space + frequency) , independent of 
N. 

For these reasons, we suspect that BEMs will lead to 
the most efficient methods to compute Casimir forces for 
complicated structures in three dimensions. Moreover, 
all that needs to be done is to take an existing BEM 
Green's function solver and change it to solve for the 



imaginary-o; Green's function. Because the imaginary-cj 
Green's functions are exponentially decaying (and ap- 
proach the familiar Poisson kernel as a; ^ 0), such a fast 
solver should actually be simpler than the correspond- 
ing real-cj solver, nor are the singularities in the Green's 
function any worse. And, as mentioned previously, the 
resulting matrix equation is real-symmetric and positive- 
definite for imaginary cj, unlike the real-o; case, even for 
dissipative materials. In short, there do not appear to be 
any substantial unsolved algorithmic problems involved 
in implementing a BEM for Casimir forces. 



IX. COMPARISON TO OTHER METHODS 

In this section, we briefly compare the stress-tensor 
approach with other known exact numerical methods ap- 
plicable to arbitrary geometries, focusing mostly on the 
computational aspects. In particular, we examine the 
methods of Emig et al. [59 and Gies et al. [58] . 

Emig's method, applicable to both separable and non- 
separable geometries and not limited to perfect metallic 
structures (although currently only demonstrated in per- 
fect metallic separable geometries), involves a surface pa- 
rameterization of the Green's function. Specifically, the 
Casimir energy is given in terms of an integral over imag- 
inary frequencies of the change in the photon density of 
states (DOS) Sp{iw): 



U 



J ^Sp{iw)dw, 



(17) 



similar to the expression in Eq. Q. However, the crucial 
aspect of this method lies in the evaluation of the DOS, 
given by [59j : 



Sp{iw) 



ld_ 

7T dw 



trln (M-iM). 



(18) 



Here M is an x dense matrix, where A' is the number 
of surface degrees of freedom, whose entries are in terms 
of the imaginary-o; vacuum Green's function G{iw] x, x') 
evaluated on the surface of each body. M^o is the same 
matrix for the case where the bodies are infinitely far 
apart. Thus, the trace in Eq. (18) is analogous to an 
integration over the surfaces of all the bodies. Inverting a 
dense matrix, multiplying two dense matrices, and taking 
the log of a dense matrix all require 0{N^) time (for 
practical algorithms) and 0{N'^) storage [90]. 

We should comment however, that Emig's method is 
closely related to a boundary-element method (BEM) as 
discussed in the previous section. BEMs also involve pa- 
rameterization in terms of surface degrees of freedom, 
which are also coupled in terms of vacuum Green's func- 
tions, leading to a dense matrix which must be inverted to 
compute the inhomogeneous Green's function. By recog- 
nizing this relationship, one should be able to exploit fast- 
multipole and similar O(A'logA') techniques to acceler- 
ate Emig's method. In particular, computations such as 
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M^M involve the solution of N linear equations, each 
of which can employ iterative methods with 0{N) stor- 
age, and similar iterative methods may also be available 
for computing matrix logarithms [90^. At best, this leads 
to O(A^^logA^) time. However, this is still much less 
efficient than the stress-tensor BEM approach discussed 
in the previous section, because the latter requires only 
^ N linear equations to be solved (the number of linear 
equations to be solved is determined by a smooth spatial 
quadrature, independent of N). 



A second exact computational method available is that 
of Ref. [SHI based on a "worldline" approach. In this 
method, the Casimir energy is represented via a scalar 
field in a smooth background potential, and the effective 
action is obtained via a Feynman path integral over all 
proper time worldlines using a Monte-Carlo approxima- 
tion. The method has only been formulated for the TM 
polarization with perfect-metal bodies, although prelim- 
inary generalizations have been suggested [91 . Specifi- 
cally, the method expresses the Casimir energy between 
two bodies as the integral of a functional 0^] over closed 
paths x(r) of length (proper time) T and center of mass 

y- 

1 rIT C 

U=-^l ^ y rfV(eE[x(r)]). (19) 

where [x(r)] is a worldline functional (similar to a step 
function) defined to be 6s[x(t)] = 1 if the path x(r) in- 
tersects a surface and otherwise [91 . This integration 
is performed by generating an ensemble of random 
A^-point paths x(r) and evaluating for each one. The 
computational complexity 0{nL ■ N - #y • #T) (where 
#y and #T are the number of y and T integrand eval- 
uations, respectively), therefore, depends on quantities 
such as the precise statistical rate of convergence of this 
integral (error ~ l/y/nZ)^ which in turn depends on the 
geometry and on the manner in which the path ensembles 
are generated and integrated. Ref. |58| does not present 
a general analysis of these quantities, nor will we do so 
here. However, in the specific case of the force between 
a radius- sphere and a plate separated by a distance 
a, Ref. 58 shows that A^ > a'^/R'^ for large a/R. This 
does not include the 0{nL ■ #y • #T) factors, although 
it seems likely that #y is at least ~ A" (so that the spa- 
tial resolutions are comparable), in which case the time 
scaling would be at least /R^ ^ N'^. In comparison, 
a BEM for the same geometry (exploiting the cylindri- 
cal symmetry to reduce it to a 2d problem similar to 
Ref. [58]) should require degrees of freedom A" that never 
scale worse than linearly with the relevant lengthscale, 
and time that scales with A" log A" rather than A'^. A 
general comparison seems difficult, however, and a de- 
tailed statistical study is outside the scope of this paper. 



X. CONCLUDING REMARKS 

The general considerations involved in designing a 
purely computational method are often quite different 
from those involved in designing an analytical method. 
For this reason, we believe it is most fruitful to start 
back at the earliest possible formulations and proceed 
using the new computational perspective, rather than at- 
tempting to add more and more corrections to analytical 
methods for specific geometries. Moreover, since decades 
of research have gone into the development of numerical 
methods for classical electromagnetism, culminating in 
methods applicable to complicated inhomogeneous three- 
dimensional geometries, it is desirable to seek approaches 
for the Casimir force that exploit these developments. 
We believe that the stress-tensor approach, developed for 
analytical calculations several decades ago, provides the 
ideal formulation for a computational approach exploit- 
ing standard numerical techniques. 

In the future, we would like to employ the stress- 
tensor approach to study Casimir forces and torques in 
more realistic and/or more unusual structures. The large 
and growing number of interesting applications of the 
Casimir effect and the ongoing experimental work on 
non-standard geometries [8 provide an environment in 
which the generality and strengths of the stress-tensor 
method could be exploited. In addition, we are currently 
implementing more efficient boundary-element versions 
of our approach in three dimensions. 

We would also like to investigate related computational 
problems. One immediate possibility is to compute the 
net torque on a body, instead of the net force. Clas- 
sically, given the stress tensor T, one can compute the 
torque by integrating r x (T • dS) [92 . This has been 
exploited by several authors to compute classical elec- 
tromagnetic torques [93, 94, 95 . Similarly, the Casimir 
torque can be obtained by using the mean (T) from 
the fluctuation-dissipation theorem, just as for the net 
force [34]. Therefore, one can compute torques by al- 
most the same method as above, via repeated evaluations 
of the classical Green's function. Several authors have 
computed Casimir torques in parallel-plate and perfect- 
metallic wedge geometries [96 , ^97\ \98\ \99\ [TOT : other in- 
teresting structures include the Casimir "pendulum" [53] 
as well as corrugated surfaces [101 , although these two 
structures have only been evaluated by methods with un- 
controlled approximations. 
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Appendix 

In what follows, we derive the force per unit length 
in translation-invariant structures in terms of an integral 
over the solutions of Bloch-periodic problems. This is not 
a new idea, but an explicit general derivation starting 
from the spatial integral of the stress tensor, including 
the case of finite periodicity, seems difficult to find in the 
literature. As given above by Eq. (13), the force per unit 
length can be written as: 



F 

L 



^ poo N-1 

-J X] JJ "^(^^5 ^ ~ ^^^) ' 



(20) 



for a periodic structure in the z-direction with period 
A, and size L = NA with periodic boundary conditions 
(ultimately, we will take the limit as ^ oo). As in 
Eq. (11), the stress tensor is expressed via the Green's 



function, and therefore Eq. (20) can be decomposed into 



individual terms of the form: 



^ POO N-1 

7 dwy2 \ \ G{i'W',rn,rn)dA, 

^ n=0 



(21) 



where = r — nAz and the Green's function is given as 
in Eq. ^ by the solution of: 



G/e(z^;x,xO = O ^J^(x-xOe/e, 



(22) 



in which O denotes the linear operator O = 
(V X V X -^w'^e). In the following, we will focus on 
the periodic z direction and leave the x and y coordi- 
nates implicit for simplicity. That is, we will write e.g. 
5{z — nA — z') instead of ^(r^ — r'). 

At this point, we can re-express the delta function over 
the periodic direction z in terms of the Fourier identity: 



6{z 



N-lN-l 



^=0 m=0 



(23) 



Substituting this into Eq. (21), we will move the 
outside the and consider the action of on the 
remaining summation: 



N-l 



Jn,m = S{z - z^ - iA - nA)e^N 



■im 



£=0 



(24) 

where we have used the periodic boundary conditions in 
L to realize that Jn,m is a cyclic shift of Jo,m with a 
phase factor. Now, we must operate on Jo,m^k and 
evaluate ai z' = z — nA. However, this corresponds to 
finding the field from a Bloch-periodic current source, 
and such a field is also Bloch-periodic. Therefore: 



(25) 

At this point, we have completely eliminated the n- 
dependence from the evaluation of the Green's function 
Jo m^ki and the phase factors from Eq. ( [24| and 
Eq. (25) cancel. The remaining summation simply 
yields which cancels the 1/A^ factor from Eq. (23). 
Equation (21) therefore becomes: 



1 rr 

- dwy2 W G{iw,m',rn,rn)dA, 

^•^0 m=0 



(26) 



where Gk{iw,m) = 0~^Jo,me/c, the field from a Bloch- 
periodic sum of delta-function sources. Finally, we can 
now take the limit A ^ oo by turning into an inte- 
gral: 



1 ^1 1 r^/A 
lim — > = / dkz 



(27) 



where kz is the Bloch wavevector (kz = 27rm/A). We 
therefore obtain Eq. (15). 
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